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The sharp change in slope of the ultrahigh energy cosmic ray (UHECR) spectrum around 10 18 ' 6 eV 
(the ankle), combined with evidence of a light but extragalactic component near and below the ankle 
and intermediate composition above, has proved exceedingly challenging to understand theoretically, 
without fine-tuning. We propose a mechanism whereby photo-disintegration of ultrahigh energy 
nuclei in the region surrounding a UHECR accelerator accounts for the observed spectrum and 
inferred composition at Earth. For suitable source conditions, the model reproduces the spectrum 
and the composition over the entire extragalactic cosmic ray energy range, i.e. above 10 17 ' 5 eV. 
Predictions for the spectrum and flavors of neutrinos resulting from this process are also presented. 


I. INTRODUCTION 

The cosmic ray spectrum spans roughly eleven decades 
of energy, 10 9 eV < E < 10 2 ° eV and has three major fea¬ 
tures: the steepening of the spectrum dubbed the “knee” 
at «sl0 15 - 6 eV [1], a pronounced hardening of the spec¬ 
trum at E & 10 18 6 eV, the so-called “ankle” feature [2- 
4], and finally a cutoff around 10 19 6 eV [3, 5]. Three ad¬ 
ditional more subtle features have been reported between 
the knee and the ankle: A hardening of the spectrum at 
around 2 x 10 16 eV [6-9] followed by two softenings at 
~ 10 16 - 9 eV [6, 7] and and 10 17 - 5 eV [2, 8-11]. The latter 
is traditionally referred to as the “second knee”. 

The variations of the spectral index reflect various as¬ 
pects of cosmic ray production, source distribution and 
propagation. The first and second knee have straightfor¬ 
ward explanations, as reflecting the maximum energy of 
Galactic magnetic confinement or acceleration capability 
of the sources, both of which grow linearly in the charge 
Z of the nucleus; the first knee being where protons drop 
out and the second knee where the highest-Z Galactic 
cosmic rays drop out. As the energy increases above the 
second knee to the ankle, the composition evolves from 
heavy to light [12] while the cosmic ray arrival directions 
are isotropic to high accuracy throughout the range [IS¬ 
IS]. Finally, as the energy increases above the ankle, 
not only does the spectrum harden significantly, but 
the composition gradually becomes heavier (interpreting 
the data using conventional extrapolations of accelerator- 
constrained particle physics models) [16, 17]. 

This observed evolution in the extragalactic cosmic ray 
composition and spectral index presents a major conun¬ 
drum. A pure proton composition might be compati¬ 
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ble with the observed spectrum of extragalactic cosmic 
rays [18] when allowance is made for experimental un¬ 
certainties in the energy scale and the fact that the real 
local source distribution is not homogeneous and contin¬ 
uous [19] (although the sharpness of the ankle is difficult 
to accommodate), but a pure proton composition is in¬ 
compatible with the depth-of-sliower-maximum (X max ) 
distributions observed by Auger [16, 17] unless current 
extrapolations of particle physics are incorrect. More¬ 
over, a fit of the spectrum with a pure proton composition 
seems to require a very strong source evolution [20] which 
leads to a predicted neutrino flux in excess of experimen¬ 
tal limits [21], On the other hand, models which fit the 
spectrum and composition at highest energies, predict a 
deep gap between the end of the Galactic cosmic rays 
and the onset of the extragalactic cosmic rays [22-27]. 
Models can be devised to fill this gap, but fine-tuning is 
required to position this new population so as to just fit 
and fill the gap [28-30]. 

Here we offer a resolution to this conundrum, by 
showing that “post-processing” of UHECRs via photo¬ 
disintegration in the environment surrounding the source, 
can naturally explain the entire spectrum and composi¬ 
tion. In our model, extragalactic cosmic rays below the 
ankle are predominantly protons from nucleons knocked 
off higher energy nuclei in the region surrounding the ac¬ 
celerator, and the spectrum and composition above the 
ankle are predominantly dictated by the accelerator and 
propagation to Earth. The model makes distinctive pre¬ 
dictions about the spectrum and flavor ratios of neutri¬ 
nos, which should enable it to be tested. If the ankle and 
the protons below it arise on account of our mechanism, 
we obtain a new constraint on UHECR sources beyond 
the Hillas criterion and total-energy-injection require¬ 
ments, namely that the environment around the source 
has the conditions giving rise to the required amount of 
photo-disintegration. 

Up until now, photo-disintegration (PD) has been 
mainly considered as a danger inside the accelerator, 
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FIG. 1. Illustration of our model calculation: Sources (yellow stars) inject cosmic rays with a power law in energy, into a 
surrounding region of radiation and turbulent magnetic fields. After propagation through this local environment and then 
intergalactic space, these cosmic rays and their spallation products are detected at Earth. The photon energies in the source 
environment are characteristically of much higher energy than in the extragalactic background light. 


as it would cut off the cosmic ray spectrum at energies 
such that the PD interaction length and the acceleration 
length are comparable. Since the acceleration length in¬ 
creases with energy, whereas the PD interaction length 
generally decreases with energy, photo-dissociation acts 
as a low-pass filter. The insight underlying the mech¬ 
anism we propose, is that if the primary locus of PD 
is outside the accelerator, PD generally acts as a high- 
pass filter, permitting the highest energy cosmic rays to 
escape unscathed while the lower energy ones are disinte¬ 
grated inside the source region, generating nucleons with 
energy 1/A of the original nucleus of mass A. As we shall 
see, these spallated nucleons naturally produce the ankle 
feature, explain why extragalactic cosmic rays below the 
ankle are protonic, and account for the spectral index 
below the ankle. Examples of systems in which the ac¬ 
celerator is embedded in a photon field and the cosmic 
rays are trapped by magnetic fields in that environment 
could be the dusty torus surrounding an active galactic 
nucleus or the interstellar medium of the star-forming 
region surrounding most young pulsars; see also [31-38]. 
The basic setup of our phenomenological model is illus¬ 
trated in Fig. 1. 

The layout of the paper is as follows. In Sec. II we 
introduce our model and in Sec. Ill we compare its pre¬ 
dictions with experimental data. Details about particle 
propagation and the calculation of multi-messenger sig¬ 
natures are given in the appendices. Section IV contains 
our conclusions. 


II. FORMATION OF THE ANKLE 

To illustrate the mechanism we have identified to cre¬ 
ate the ankle and generate protons below, consider a 
system in which the accelerator (also referred to as the 
source) is embedded in an environment in which the cos¬ 
mic rays are confined for some time by magnetic fields 
while interacting with the ambient radiation held. Our 
essential simplifications are: (i) a fast acceleration mech¬ 


anism and/or a low photon density inside the accelerator, 
(ii) no energy is lost except through an interaction, and 
whenever a nucleus interacts it loses one or more nucleons 
by photo-disintegration or photo-pion production (in this 
case the nucleus loses a fraction of its energy correspond¬ 
ing to the reduction in its nuclear mass); (Hi) a cosmic 
ray either escapes without changing energy, with a rate 
T esc , or the cosmic ray interacts one or more times before 
escaping; (iv) r esc and T- lnt are independent of position 
in the source environment and depend only on {E, A, Z} 
of the nucleus. In this approximation the number of nu¬ 
clei in a given energy range and with a specified {A, Z} 
decreases exponentially with time, with 

t = + rr ])- 1 . ( 1 ) 

A fraction 

lesc = (1 T T eS c/T itl t) (2) 

of the particles escape without interaction and the rest 
interact before escaping, so r/i nt = 1 — i] esc . Note that 
Tjesc and 77i nt depend only on the ratio of the escape and 
interaction times, but not on the absolute value of either 
of them. 

A simple analytic treatment is instructive. To illus¬ 
trate the low/high-pass filter mechanism, consider the 
case that the escape and interaction times are both power 
laws in energy, 

T esc = a {E/E 0 ) s and r int = b (E/E 0 ) c . (3) 

Then 

Vesc(E) = (1 + R 0 (E/Eo) 5 ^)- 1 , ( 4 ) 

where Rq = a/b is the ratio of the escape and inter¬ 
action time at reference energy Eq. When S > </, the 
source environment acts as a low-pass filter on the parti¬ 
cles injected from the accelerator, leading to a cutoff in 
the escaping spectrum at high energies. This situation 
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FIG. 2. Interaction times of 28 Si in a broken power-law pho¬ 
ton field with parameters a = §, /3 = — 1 and £o = 0.11 eV. 
Top panel: photo-disintegration, middle panel: photo-pion 
production, bottom panel: sum of the two processes. The 
results of numerical integration using detailed cross sections 
are shown as thick solid lines, while those of the narrow- 
resonance-approximation (detailed in Appendix B) are dis¬ 
played with thin dashed lines. 


is typical of leaky box models of diffuse acceleration at 
time-independent shocks [39-41] where <5 > 0 because the 
higher the energy of the particle, the longer it needs to 
stay in the accelerator to reach its energy. By contrast, 
if the escape time decreases with energy, as in the case 
of diffusion in turbulent magnetic fields outside the ac¬ 
celerator, then it is possible to have S < £ leading to a 
high-pass filter on the energy spectrum of injected nuclei: 
the lower the energy, the more time the nuclei have to 
interact before escaping, leading to a hardening of the 
spectrum and lightening of the composition of nuclei es¬ 
caping the region surrounding the source. The spallated 
nucleons have energies of E = E^/A-, these nucleons are 
most abundant at low energies and have a steeper spec¬ 
trum oc (1 — rj esc (E* A 1 )). Thus the high-pass scenario 
leads naturally to an ankle-like feature separating the 
nucleonic fragments from the remaining nuclei. The nor¬ 
malization and slope of the spectrum of spallated nucle¬ 
ons relative to that of the primary nuclei is determined 
by how thoroughly the primary nuclei are disintegrated, 
which is governed by the ratio of escape and interaction 
lengths of the most abundant primaries. 


To obtain a more realistic treatment of the interac¬ 
tion time, we must specify the shape of the spectrum of 
the target photons. In our work to date we have con¬ 
sidered: (i) a broken power-law (BPL), characterized by 
its peak energy eo and lower and upper spectral indices 
a, (3 (this is a simplified representative of non-thermal 
emission that allows for analytic calculation as discussed 
below and in Appendix B); (ii) a black-body spectrum; 
(in) two types of modified black-body spectrum, which 
result from a reprocessed black-body in a dusty environ¬ 
ment [42]. Details are given in Appendix A. For such 
peaky photon spectra the interaction time does not have 
the simple representation of (3) but it does have a rather 
universal structure. In our actual calculations we adopt a 
numerical integration of Talys [43, 44] and Sophia [46] 
cross sections using [47], but the analytic expression for 
T; nt derived in Appendix B for the BPL in the narrow- 
width approximation for the interaction cross sections, 
is qualitatively similar and useful for understanding. As 
can be seen in Fig. 2 the folding of a single resonance with 
a broken power-law spectrum leads to a “V” shape curve 
for Tint hr a log-log plot for both photo-disintegration (top 
panel) and photopion production (middle panel). Com¬ 
bining both processes in narrow-resonance approxima¬ 
tion yields an interaction time with a “W” shape, while 
numerical integration including the plateau for nrulti- 
pion production softens the “W” to what we shall refer to 
as an “L” shape for brevity, a shown in the bottom panel 
of Fig. 2. As evident from Fig. 2, below the inflection 
point for photodisintegration E b , the narrow-resonance 
approximation is good, while from the full numerical in¬ 
tegration in the high-energy region Tint is roughly con¬ 
stant, so using the BPL spectrum, we have the approxi¬ 
mate representation: 


Ti„t(-E) « n 


(E/E b y +1 E^E b 
1 E > E b ’ 


( 5 ) 


where formulae for r b and E b are given in Appendix B, 
and the parameter values for photodisintegration are to 
be used. 

Returning to the discussion of Ti nt in (3) with (5) yields 
the fraction of nuclei which escape without interaction in 
a peaky photon spectrum. It is straightforward to see 
that if 6 < 0 and the interaction time is described by an 
L-shaped curve, then r) esc has the properties of a high- 
pass filter. These conclusions do not depend on the exact 
shape of the photon spectrum. As can be seen in Fig. 13 
of Appendix A, the interaction times flatten to an L- 
curve as well if the photon density is assumed to follow 
a (modified ) black body spectrum. 


III. COMPARISON WITH EXPERIMENT 
A. Fiducial Model 

As our fiducial example, we adopt a broken power-law 
photon spectrum as a simplified representative of non- 
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thermal emission given by 


n(e) 


= n 


BPL 

0 


(e/eo)“ £ <£o 
(e/eo ) /3 otherwise. 


( 6 ) 


where e is the photon energy, the maximum photon num¬ 
ber density is at an energy of e 0 and following [39] we take 
the slope parameters a = + § and /3 = — 2 . As we shall 
see later, any peaky spectrum gives similar results, with 
the position of the peak, eq, being the most important 
parameter besides the peak photon density. 

Inspired by the energy dependence of the diffusion co¬ 
efficient for propagation in a turbulent magnetic field, we 
model T esc as a power law in rigidity E/Z, 

T eS c = MEZ-'/EoY. (7) 

Since only the ratio of escape and interaction times mat¬ 
ters, and the {E, A , Z } dependence of this ratio is entirely 
determined once the spectral index of the escape time <5 
is specified, the remaining freedom in characterizing the 
source environment can be encoded by specifying the ra¬ 
tio of escape to interaction time for a particular choice 
of {E, A, Z} : which we take at 10 19 eV for iron nuclei, 
denoted In application to a particular source can¬ 
didate, Rig depends on the density of photons and the 
properties of the turbulent magnetic field that delays the 
escape of the UHECRs from the environment of their 
source. 

Figure 3, upper panel, shows the escape and interac¬ 
tion times in the fiducial source environment, as a func¬ 
tion of the cosmic ray energy, for proton, He, N, Si and 
Fe; the interaction times are calculated including both 
photo-disintegration and photo-pion production. The 
gross features of the energy dependence of the interaction 
times can be understood in the approximation of reso¬ 
nant interactions in the nucleus rest frame £' es . At low 
cosmic-ray energies, reaching £' les requires high photon 
energy (e > eq), so that the interaction time decreases 
with increasing cosmic-ray energy as r oc E^ +1 . However 
for high enough cosmic ray energy, the resonance can be 
reached in collisions with photons of e < e 0 . From here, 
as the cosmic ray energy increases, the photon density at 
the resonant energy decreases as e“, and correspondingly 
the interaction times increase. The laboratory energy of 
the inflection point of the interaction times for a cosmic 
ray nucleus of mass Am p is at E = Am p £[. es /(2£ 0 ). The 
inflection point of the photo-dissociation times can be 
seen as a dip in the plot in the upper panel of Fig. 3, e.g., 
at around 10 188 eV for iron nuclei. At slightly higher 
energy, photo-pion production becomes important, with 
the result that the energy dependence of the interaction 
time is roughly speaking an L-shaped curve in a log-log 
presentation. 

Using these energy-dependent interaction and escape 
times, we propagate nuclei through the source environ¬ 
ment with the procedure described in Appendix C. Cos¬ 
mic rays of some given composition are injected from the 
accelerator into the source environment with a power law 
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FIG. 3. Top: Interaction and escape times for A = 
1, 4, 14, 28 and 56 (bottom to top for escape; vice versa for 
interaction) for the fiducial model photon field with eo = 0.11 
eV. Bottom: Injected 28 Si flux (bold dashed) and escaping 
fluxes: thin black solid line denotes the sum of all escap¬ 
ing nuclei and solid curves give contribution of different mass 
groups with low energy intercept increasing with mass. Nu¬ 
cleons from photo-dissociation and photo-pion-production are 
shown with thin-dashed and dotted curves, respectively. 


spectrum and an exponential cutoff at some maximum 
rigidity. To keep the complexity of the fiducial model to 
a minimum, we inject only a single nuclear species and fix 
the injection spectral index 7 = — 1 , as expected for accel¬ 
eration in young neutron stars [48]. The particles escap¬ 
ing the source environment are then propagated through 
the intergalactic medium using the procedure explained 
in Appendix D. 

In total the fiducial model has 14 parameters, with 8 
parameters allowed to float freely in the fit, as indicated 
in Table I. The spectral index and normalization of the 
Galactic spectrum are free “nuisance” parameters with 
the best fit giving a spectral index of —4.2. This should 
be understood as an effective spectral index describing 
the cutoff of the Galactic cosmic ray population, and 
hence cannot be directly compared with the parameter 
reported by the KASCADE-Grande Collaboration [49], 
because their single-power law fit is driven by the “low- 
energy” data. The fraction of Galactic cosmic rays at 
10 17 ' 5 eV is 55%. 

The best description of the data is obtained with 28 Si 
of maximum energy Z 10 18 5 eV = 4.6 x 10 19 eV; the im¬ 
pact of allowing other parameters to vary is discussed 
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in following sections. Normalizing this model to the ob¬ 
served flux at Earth, we infer a comoving volumetric en¬ 
ergy injection rate in CRs at z = 0, above 10 175 eV, of 
ei 7.5 = 9.2 x 10 44 erg Mpc -3 yr -1 . 

The unmodified injection spectrum and the spectrum 
of escaping nuclei for this fiducial model are shown in 
the lower panel of Fig. 3. At low energies, the nuclei are 
depleted relative to the injected flux because T esc » Tint, 
but the escaping nuclei follow the original spectral index 
because in this example the interaction and escape times 
are parallel, as to be expected for 5 = /3 + 1. Once the 
corner of the L-shape is reached, the fraction of escaping 
nuclei grows, leading to an apparent hardening of the 
spectral index. 

Even for the simple case in which a single nuclear 
species is injected into the source environment, we ob¬ 
tain a complex evolution of the mass composition with 
energy. At low energies the composition is dominated by 
knock-off nucleons whereas at high energies the composi¬ 
tion becomes heavier as the ratio of escape to interaction 
time drops and more heavy nuclei can escape before in¬ 
teracting. 

This fiducial model of interactions in the source envi¬ 
ronment is a very simple one, yet even so it offers a re¬ 
markably good accounting for the flux and composition 
at Earth as determined by the Pierre Auger Observatory. 
(Data from the Telescope Array (TA) are consistent with 
the Auger results within systematic and statistical un¬ 
certainties [51, 52] and also can be well-fit; we come to 
TA separately below.) In Fig. 4 we compare the fidu¬ 
cial model prediction to the Auger measured flux, from 
10 17 ' 5 eV to above 10 20 eV [50] and to the mean and vari¬ 
ance of the distribution of the logarithm of mass on top 
of the atmosphere, (In A) and V(ln A) [16, 53, 54]. There 
is a good overall agreement between the model and the 
data. The shape of the spectrum is described well, in¬ 
cluding the ankle and the flux suppression. The model 
also qualitatively reproduces the increase of the average 
logarithmic mass with energy and the decrease of its vari¬ 
ance. 

The neutrino signals of the fiducial model are shown 
in Fig. 5; details of the calculation are given in Appendix 
E. An exciting aspect of our model for the ankle is the 
presence of a detectable anti-electron-neutrino flux from 
neutron /3-decay, with a rate consistent with the naive es¬ 
timate of [57] . The right panel of Fig. 5 shows the number 
of events as a function of energy predicted in ten years of 
IC86, using the IceCube acceptance for different neutrino 
flavors given in [56]. In total, the fiducial model predicts 
3.5 events in the range 10 16 — 10 1 ' eV after 10 years of 
operation of IceCube (corresponding to about one year of 
operation for an upgraded IceCube-Gen2 detector [58]). 
We emphasize the distinctive v e enrichment due to beta 
decay of spallated neutrons. 

The associated photon flux from nuclear de-excitation 
in our model is well below the Fermi-LAT data (see 
Appendix E for more detail). Photo-pion interactions 
at the source and during propagation produce an addi¬ 


tional flux of photons via 7r°-decay; this is consistent with 
Fermi-LAT data, as follows: If the origin of the photons 
measured by Fermi-LAT is exclusively from these inter¬ 
actions, then from [59] the associated diffuse neutrino 
flux saturates the IceCube upper limit [56]. Since the 
neutrino flux in the fiducial model is below the IceCube 
limit, it follows that also the associated photon flux is 
consistent with Fermi-LAT data. A more sophisticated 
realization of our mechanism than in the fiducial model 
must also respect the IceCube limits, and therefore the 
Fermi-LAT data as well. 


B. Model Variations 

In this section we discuss the impact of theoretical and 
experimental uncertainties on our model, as well as dif¬ 
ferent choices for the fiducial parameters. 

1. Experimental Uncertainties 

To study the influence of the experimental systematic 
uncertainties on our fit, we have repeated the fit for all 
combinations of altering the measurements by +1, +0 
and — lcr sys . of the quoted uncertainties on the energy 
and composition scale. We find that the best fit is ob¬ 
tained within the experimental systematics when shifting 
the energy scale up by + lcr sys . = +15% and by shifting 
(In A) and V (In A) corresponding to a shift of the shower 
maximum by — l<j sys . —10 g/cm . The best-fit values 
after the application of these shifts are shown in brackets 
in Table I. Most notably, the peak energy of the photon 
spectrum decreases from 110 to 70 meV and the best-fit 
value of the spectral index of the escape time decreases 
from ~ —3/4 to almost —1. The neutrino flux at Earth 
obtained for this fit is about 30% smaller than in case of 
the fiducial model. This is mainly due to the difference in 
the best-fit peak energy of the photon field in the source 
environment. The sensitivity of the neutrino flux to £q 
will be further discussed in Sec. Ill B 5. 

The overall description of the spectrum and composi¬ 
tion is considerably improved, as can be seen in Fig. 6. 
The model variations discussed below will therefore be 
performed based on shifted data. 

2. Hadronic Interactions in Air Showers 

The interpretation of experimental air shower data in 
terms of mass composition relies on the validity of extrap¬ 
olations of the properties of hadronic interactions to ul- 
trahigh energies. Using alternative models for this inter¬ 
pretation (Sibyll2.1 [62] or QGSJetII- 04 [63] instead 
of Epos-LHC [64]), decreases the value of the (In A) data 
points by about (In A) = —0.6 and leads to a worse fit 
of the data. If this difference between models gives a fair 
estimate of the uncertainties of the mass determination 
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(b) Composition at Earth 


FIG. 4. Spectrum and composition at Earth. The data points are from the Pierre Auger Observatory [16, 50], error bars denote 
the statistical uncertainties and the shaded boxes illustrate the experimental systematic uncertainties of the composition. The 
composition estimates are based on an interpretation of air shower data with Epos-LHC; the lines denote the predictions of 
our fiducial model. 
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FIG. 5. Neutrino spectrum (left) and expected number of events in 10 IC86-years (right) for the fiducial model. The measured 
flux of low-energy extragalactic neutrinos from IceCube [55] is shown in the left panel (purple lines) as well as the 90% CL 
upper limit on the flux of high-energy neutrinos (dashed area) [56]. The peak in the electron neutrino flux at about 10 1;> ' 8 eV 
seen in the right panel is due to the increased interaction probability of anti-electron neutrinos at the Glashow resonance. 


in both directions, er t heo((ln4)) = +0.6, then a hadronic 
interaction model that leads to a heavier interpretation 
of Auger data than EPOS-LHC would make the fit with 
the fiducial model even better, similar to the systematic 
shift in the composition scale discussed in the previous 
section. 


3. Mass Composition at the Source 

It is remarkable that a good description of both the 
spectrum and mass composition at Earth is possible by 


assuming only a single injected species at the source as 
assumed for simplicity in the fiducial model. However, 
depending on the astrophysical scenario, this might be 
an unrealistic assumption. 

In Fig. 7 we explore the capability of our model to in¬ 
corporate additional flux components of mass A\ below 
and above the mass A 2 ~ 29 that gives the best fit for the 
fiducial single-mass model. As can be seen, our calcula¬ 
tion allows for an additional proton or helium component 
as large as 80% and up to 70% for nitrogen. 

For an additional flux component with a heavy mass, 
the model is more restrictive as illustrated in the lower 
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source parameters 

power law index of injected nuclei 

7 

fix 

-1 

mass number of injected nuclei 

A 

free 

28 (29) 

maximum energy 

E p 

^max 

free 

10 !8.5 (18.6) e y 

cosmic ray power density, E > 10 17 ' 5 eV 

£17.5 

free 

9.2 (13) xlO 44 erg Mpc -3 yr~ 4 

evolution 

£(z(t)) fix 

star formation rate [60] 

source environment 

energy of maximum of photon field density 

£0 

free 

0.11 (0.07) eV 

power law index of photon spectrum (e < £ 0 ) 

a 

fix 

+ 3 
' 2 

power law index of photon spectrum (e ^ £ 0 ) 

p 

fix 

-2 

power law index of escape length 

s 

free 

-0.77 (-0.94) 

ratio of interaction and escape time 

TD Fe 
-^19 

free 

4.4 (3.7) xlO 2 

propagation to Earth 

infra-red photon background 

- 

fix 

Gilmorel2 [61] 

spectrum of Galactic cosmic rays 

power law index at Earth 

7gal 

free 

-4.2 (-3.7) 

mass number of Galactic nuclei 

Agal 

fix 

56 

flux fraction at 10 17 ’ 5 eV 

/gal 

free 

57 (72) % 


TABLE I. Parameters of the fiducial model. Values in brackets denote the parameters of the best-fit obtained when shifting 
the data by its systematic uncertainties (see text). 


left panel of Fig. 7 using A± = 56. In this case, the de¬ 
scription of the data considerably deteriorates for frac¬ 
tions above 10%. The reason for this behavior is twofold. 
Firstly, the injection of too much iron at the source leads 
to a too heavy composition at Earth as compared to the 
estimates from the Pierre Auger Observatory. Secondly, 
if the end of the cosmic ray spectrum is to be described 
by the maximum rigidity of iron nuclei, then the energy 
of secondary nucleons needed to populate the flux at and 
below the ankle is too small to describe the data (the 
maximum energy of secondary nucleons is 1 /A of the 
maximum energy of nuclei). 


If the cut-off of the flux is at higher energies, as sug¬ 
gested by the measurement of TA [65] , then a larger frac¬ 
tion of iron primaries at the source can be incorporated, 
as shown in the lower right panel of Fig. 7. When using 
the TA data in the fit, as shown in Fig. 8 , the spectrum 
can be described reasonably well even for an injected 
flux consisting of 100% iron nuclei. But in this scenario 
the composition at Earth at ultrahigh energies is heavier 
than suggested by the interpretation of the X max data of 
Auger. 


As an illustration of a more complex composition 
model, we use the abundances of Galactic nuclei at a 
nucleus energy of 1 TeV, which we read from Fig. 28.1 
in [ 66 ]. The flux fractions are 0.365, 0.309, 0.044, 0.077, 
0.019, 0.039, 0.039, 0.0096, 0.014, 0.084 for H, He, C, 
0,Ne, Mg, Si, S, Ar+Ca, Fe, respectively. The result¬ 
ing fit is shown in Fig. 9 (7 = —1.25 and 5 = —1). This 
example demonstrates that our mechanism for producing 
the ankle is working even when considering a complicated 
mix of primaries. 


4- Source Evolution and Spectral Index 

To have a concrete fiducial model, we needed to specify 
how the production of UHECRs varied over cosmological 
time scales. This is known as the source evolution, which 
we took to be in direct proportion to the star-formation- 
rate - as would be expected in a source scenario such as 
young magnetars. In this section, we consider alternative 
evolutions of the source luminosity density described by 
the simple one-parameter functional form 




(1 + z) m Z< Zo , g , 

(1 + z 0 ) m exp (—(2 — 20 )) otherwise 


with zo = 2 and m ranging from —4 to +4. m = 0 
yields a uniform source luminosity distribution, m = +4 
corresponds to a strong evolution similar to the one of ac¬ 
tive galactic nuclei, and negative values result in sources 
that are most abundant or most luminous within the low- 
redshift universe as suggested in [67]. The resulting fit 
parameters are displayed in Fig. 10 for three choices of 
the spectral index 7 of the injected flux: — 1 , as in the 
fiducial model, —2 as expected for stochastic shock accel¬ 
eration and for letting 7 float freely in the fit. As can be 
seen in Fig. 10(a), 7 = —2 gives a poor description of the 
data for in > 0 , but is a viable choice for closeby sources, 
in accordance to the findings of [67]. For positive values 
of m, a fixed value of 7 = —1 gives a similar fit quality 
as the freely floating 7 , but the latter converges to val¬ 
ues larger than —1 for source evolutions with m > 2 (cf. 
Fig. 10(c)). 

For the “traditional” source evolutions with to 3 s 0 and 
the fit with 7 = — 1 we find that most of the parameters 
exhibit only a minor variation with to, with the exception 
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(a) Tj t | i and Tesc ■ 



(b) Injected (dashed line) and escaping (solid lines) fluxes. 





lg(E/eV) lg(E/eV) lg(E/eV) 

(c) Flux at Earth (d) Composition at Earth 


FIG. 6 . Spectrum and composition at Earth. The data points are from the Pierre Auger Observatory [16, 50] shifted, by plus 
one sigma of systematic uncertainty for the energy scale and minus one sigma for the X max scale. The lines denote the best-fit 
within our fiducial model. 


of the power-law index of the escape time S (Fig. 10(e)) 
and the power density £ 17.5 (Fig. 10(e)). 

We conclude that our model for the ankle does not 
critically depend on the choice of the source evolution, 
but that for a given choice of m we can constrain the 
allowed values of 7 , 6 and £ 17 . 5 . 


5. Photon Spectrum 

We repeated the model fits using alternative energy 
distributions of the photon density instead of the broken 
power law used in the fiducial model: a black body spec¬ 
trum and two modified black body spectra. All four spec¬ 
tra are normalized to the same integral photon density 
and depend only on one parameter, the peak energy £0 
(see Appendix A). The resulting fit results are shown in 
Fig. 11 for a freely floating spectral index 7 and for source 
evolutions with m O 0. As can be seen, all four photon 
spectra describe the data equally well (Fig. 11(a)). The 


best-fit values of the free model parameters are very sim¬ 
ilar and in particular the obtained peak values are within 
+20 meV. We conclude that as long as the photon spec¬ 
trum is “peaky”, the particular details of its shape do 
not influence the parameters of our model. 

The sensitivity of the fit to the peak energy is shown in 
the left panel of Fig. 12. As can be seen, the x' 2 deterio¬ 
rates very quickly at low values of £ 0 , but it is almost flat 
above the minimum. This feature can be easily under¬ 
stood recalling e 0 in the “L-curve” approximation intro¬ 
duced in Sec. II: The smaller e 0 , the larger is the energy 
of inflection of the interaction length, £),. For too-small 
values of £o; the interaction and escape times are parallel 
over the full energy range and thus no high-pass filter is 
created. On the other hand, once Eb is small enough, 
a further decrease changes only the flux at low energy, 
where the escaping spectrum is dominated by low-mass 
nuclei from spallation (see e.g. Fig. 3) which can be com¬ 
pensated by adjusting other parameters such as 
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(a) Proton. 


(b) Helium. 


(c) Nitrogen. 



(d) Iron, Auger flux. 


(e) Iron, TA flux. 


FIG. 7. Injection of two mass components. The first mass value, Ai, is fixed and contributes the fraction indicated on the 
x-axis to the total flux. The second mass value, A 2 , is varied as shown on the y-axis. The fit quality is indicated by the colors. 


To first order, our model can therefore only give a lower 
limit on the peak energy of the photon flux in the source 
environment. However, future limits or observations of 
neutrinos in the 10-100 PeV range will help to constrain 
this important source property, because the number of 
predicted neutrinos strongly depends on £ 0 , as shown in 
the left panel of Fig. 12 by the superimposed open sym¬ 
bols. A larger peak energy of the ambient photon en¬ 
vironment increases neutrino production at the source 
in two ways. Firstly, shifting Ei, to lower energies (and 
compensating as necessary by adjustment of Rig) moves 
the interaction times of protons closer to the escape time 
and correspondingly additional neutrinos are produced 
via photo-pion production of protons (compare e.g. the 
red curves at around 10 18 eV in the upper panel of Fig. 3 
(£0 = 110 meV) to the ones in Fig. 6 (a) (£0 = 70 meV)). 
Secondly, increasing e 0 moves the minimum of the inter¬ 
action time for photo-pion production of nuclei to lower 
energies. Since the neutrinos from photo-pion production 
carry a larger fraction of the nucleon energy than the neu¬ 
trinos from neutron decay after photo-dissociation, this 
increases the neutrino flux as well. 


It is tempting to give a quantitative interpretation 
of the x 2_cul ' ve °f Fig - 12 in terms of a lower limit on 
£0 and the number of neutrinos. However, the mini¬ 
mum of x 2 is far away from x 2 /-^df = 1 which - as¬ 
suming this model is correct - is indicative of experi¬ 
mental systematics or an under-estimation of the experi¬ 
mental uncertainties or of deficiencies in the modeling of 
hadronic interactions in the atmosphere needed to inter¬ 
pret the data in terms of mass composition (see above). 
In the absence of a concrete explanation we follow the 
PDG [ 66 , 68 ] and rescale the uncertainties by a common 
factor S = (Xmm/^df) 5 1° bring the rescaled x 2 /-^df 
to 1. This rescales the y 2 value of any given model so 
that the number of standard-deviations it is from the 
minimum is given by N' a = S'* 1 Vx'modei ~ Xmin- This 
yields an approximate lower limit on £ 0 at N' a = 3 of 
£0 > 34 meV and N v ( 10 x IC 86 ) > 0.4 assuming the va¬ 
lidity of the fixed fiducial parameters given in Table I. 
The corresponding lower temperature limits are 180 K, 
125 K and 100 K for the black body spectra with a = 0, 
1 and 2 respectively. The lower limit on the neutrino 
spectrum is shown in the right panel of Fig. 12. 
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(a) Tj n t and. TeSC • 
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(b) Injected (dashed line) and escaping (solid lines) fluxes. 



lg(E/eV) 



(c) Flux at Earth 


(d) Composition at Earth 


FIG. 8. Spectrum and composition at Earth. The data points are from the TA [65] (flux) and the Pierre Auger Observatory [16] 
(composition). The latter have been shifted in energy to match the energy scale of TA and the A' max scale is shifted down by 
1 sigma. The lines denote the fit with our model assuming a pure iron composition at the source. 


6. Hadronic Interactions in the Source Environment IV. CONCLUSIONS 


In addition to interactions with the background pho¬ 
ton field, nucleons and nuclei can also scatter off hadrons 
in the source environment. In this paper we assume 
that the density of hadronic matter in the source en¬ 
vironment is low enough that such hadronic interactions 
can be neglected. For any concrete astrophysical real¬ 
ization of our scenario, one must check and if necessary 
include hadronic interactions in the source environment. 
Production of tt^’s and 7r°’s in hadronic collisions could 
significantly increase the fluxes of neutrinos and photons 
emitted in the EeV energy range. Fast-spinning newborn 
neutron stars provide a particular example [69]. Pre¬ 
cise estimates of the impact of hadronic collisions on the 
predictions of our model will be presented in a separate 
publication. The results presented here are valid for all 
astrophysical systems in which the interactions are dom¬ 
inated by photo-nuclear processes. 


In this paper we have proposed a new explanation for 
the ankle in the cosmic ray spectrum, and for the evo¬ 
lution with energy of the composition of extragalactic 
cosmic rays: from light below the ankle to increasingly 
heavy above. When nuclei are trapped in the turbulent 
magnetic field of the source environment, their escape 
time can decrease faster with increasing energy than does 
their interaction time. Under these conditions, only the 
highest energy particles can escape the source environ¬ 
ment unscathed, and the source environment acts as a 
high-pass filter on UHECRs. Nuclei below the crossover 
energy such that T esc > T; n t interact with photons in the 
environment around the source, with ejection of nucle¬ 
ons or alpha particles and consequent production of a 
steep spectrum of secondary nucleons. The superposition 
of this steeply falling nucleon spectrum with the harder 
spectrum of the surviving nuclear fragments creates an 
ankle-like feature in the total source emission spectrum. 
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(a) Injected fluxes in 5 mass groups. 


(b) Escaping fluxes (sum of injection shown as dashed line). 




(d) Composition at Earth 


FIG. 9. Spectrum and composition at Earth. The data points are from the Pierre Auger Observatory [16, 50] shifted by their 
systematic uncertainty as in Fig. 6. The injected composition follows a Galactic mixture with 10 elements (see text). 


Above the ankle, the spectrum emerging from the source 
environment exhibits a progressive transition to heavier 
nuclei, as the escape of non-interacting nuclei becomes 
efficient. Abundant production of D e 's is a signature of 
this mechanism. 

We illustrated the high quality of the fit which can 
be obtained to the Auger data, with a fiducial model 
in which nuclei are accelerated up to a maximum rigid¬ 
ity found to be as 10 185 V, with spectrum ccE~ x , and 
are then subject to photo-disintegration in the vicinity 
of the accelerator before escaping for their journey to 
Earth. We showed that the details of the photon spec¬ 
trum around the accelerator are unimportant, except for 
its peak energy. The other important characteristic of 
the environment is the photon density relative to the 
magnetic diffusivity, which we characterized in a very 
simplistic way (through a single parameter) in this ini¬ 
tial study. We studied the sensitivity of the mechanism 
to the energy-scale uncertainty and hadronic-interaction¬ 
modeling uncertainty, which affects the composition in¬ 


ferred from the atmospheric shower observations, and 
also used the TA spectrum instead of the Auger spec¬ 
trum. The conclusion of these studies is that a good 
quality fit can be obtained in most cases, but details of 
the fit parameters such as the composition and maximum 
energy characterizing the accelerator change. A corollary 
is that until these systematic uncertainties in the obser¬ 
vations and their interpretation are reduced, such details 
of the accelerator cannot be reliably inferred from the 
data. The fiducial model parameters needed in the fits 
are such that the scenario can be reasonably achieved in 
at least one type of proposed astrophysical source, as will 
be discussed in a future publication. 

Our mechanism has two predictions beyond fitting the 
shape of the spectrum and composition evolution, which 
are independent of many environmental variables and can 
be used to test the validity of this scenario for produc¬ 
tion of the ankle, i) The spectral cutoff of spallated nu¬ 
cleons emerging from the source environment is 
where i? max is the rigidity cutoff of the accelerator, be- 
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(c) spectral index of injected spectrum. (d) ratio of escape and interaction time. 
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(g) injected mass 


(h) UHECR energy injection rate 


FIG. 10. Fit results as a function of source evolution for different spectral indices of the injected flux: 7 fixed to —1 (open 
squares), fixed to —2 (open circles) and best fit (filled circles). On the x-axis the power m of the source evolution is shown; the 
last bin reports the values for the fiducial model (SFR) evolution from [60], Eq. (D 6 ). 
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FIG. 11. Fit results as a function of source evolution for different photon spectra: Broken power law (BPL, open squares), 
black body spectrum (BB, open circles), modified black body spectrum (MBB) with a = 1 (filled circles) and a = 2 (filled 
squares). On the x-axis the power m of the source evolution is shown and in the last bin the fit values for the fiducial evolution 
from [60], Eq. (D6), is shown. 















14 



£ o [meV] 


— v e -- - v, —sum 53IC2013 

— v e — v„ — v, — IC2014 



lg(E/eV) 


FIG. 12. Left: Fit quality of the fiducial model (closed symbols) and number of neutrinos (open symbols) as a function of 
peak energy eo of the photon spectrum in the source environment. Four types of photon spectra are shown: Broken power law 
(BPL), black body spectrum (BB) and two modified black body spectra (MBB). The minimum y 2 °f BPL corresponds to the 
result shown in Fig. 6. Right: Lower limit on the neutrino flux obtained for a modified black body spectrum with a = 2, and 
£o = 34 meV (T = 100 K). The lines and hatched area at the top of the figure are the measured neutrino flux and upper limit 
from IceCube [55, 56] 


cause f?max, spal .nuc. -E'ma x,a/ -A while -E'max,A — Z R max? 
and finally Z/A = largely independent of composition. 
This relation holds prior to the extragalactic propaga¬ 
tion from the source, thus giving complementary infor¬ 
mation on the accelerator to that obtained from the spec¬ 
trum and composition above the ankle alone, ii) There 
is a one-to-one relation between the spectrum of spal- 
lated nucleons and the anti-electron-neutrinos produced 
by beta decay of neutrons, unless the spallated nucle¬ 
ons lose energy by interacting with hadronic material in 
the source environment. Independent of other proper¬ 
ties of the environment or the source evolution, P e ’s will 
have an identical spectral shape, shifted down by a factor 
~ 1/1000 from the kinematics of n —> pe~v e and reduced 
by a factor-2 in normalization because only half the nu¬ 
cleons are neutrons. This follows because propagation 
energy losses are small for nucleons of such low energy, 
and redshift impacts both nucleons and neutrinos identi¬ 
cally. Thus, detailed comparison of the v e and spallated 
nucleon spectra will reveal if hadronic interactions in the 
source environment are important, which would imply a 
correlated production of photo-pion produced neutrinos. 


exploring another mechanism for producing the ankle, 
arising in the context of gamma-ray bursts [70]. 
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After this work was presented at the IceCube Particle 
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FIG. 13. Left: Comparison of photon spectra. BPL: Broken power law (solid), BB: black body spectrum (dashed), MBB: 
modified black body spectrum (dotted and dash-dotted). The curves are normalized to match the integral of the black body 
spectrum and the temperatures are chosen to match the peak energy of the broken power law. Right: Interaction times 
corresponding to the four photon spectra. 


Appendix A: Photon spectra 


In this paper we explore the propagation effects of the four types of photon spectra shown in Fig. 13. The first 
consists of broken power-law (see e.g. [39]) as a simplified representative of non-thermal emission given by 


n(e) 


= n 


BPL 

0 


(e/e 0 )“ £ < £ 0 
(e/eo)^ otherwise. 


(Al) 


where e is the photon energy and the maximum of the number density is at an energy of Eq. 
We also consider modified black-body spectra using the functional form 


i{e) = 


MBB 


8tt 


(he) 3 e^T — 1 


(A2) 


where T denotes the temperature in the case of pure black-body, and the absorption factor is given by (j^) cr (see 
e.g. [42]). h, k and c are the Planck constant, Boltzmann constant and speed of light respectively. For a = 0 and 
77 MBB _ 1 unmodified black-body spectrum is obtained. The relation between the peak energy and temperature 
parameter is given by a modified Wien’s displacement law, 

£ 0 = [w(-e- & b ) + b] k T, (A3) 

where W(x) is the Lambert function (see e.g. [71]) and b = a + 2. 


For the study of the effect of using different functional forms of photon spectra in our model (cf. Sec. Ill B 5), it is 
useful to use a common normalization for all spectra. The integral photon density of Eq. (Al) is 


Ibpl = n-Q PL £q 


\a + 1 



and for Eq. (A2) it is 


Z MBB 


= n; 


MBB 


87 r 
(he) 3 


(fcT) 3 C(o- + 3,l)r(o- + 3), 


(A4) 


(A5) 


where £(x) denotes the Riemann zeta function and F(x) is the Gamma function. Choosing the photon density of the 
unmodified black body spectrum as reference we use the following normalization constants, 

n o PL = -^mbb/^bpl 


(A 6 ) 
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and 

"TV) = ^mbbAmbb- (A7) 

An example of the four photon spectra after normalization and for the same peak energy of £q = 50 meV is shown 
in the left panel of Fig. 13. The corresponding interaction time for the sum of photo-dissociation and photo-pion 
production is shown in the right panel. 


Appendix B: Photo-Nuclear Interactions 


The interaction between photons and high energy nuclei has been extensively discussed in the literature [72-82] . In 
this appendix, we describe how we implement the photon-nucleus collisions in our analysis. The interaction time for 
a highly relativistic nucleus with energy E = 'yAirip (where 7 is the Lorentz factor) propagating through an isotropic 
photon background with energy e and spectrum n(e), normalized so that the total number density of photons is 
\n(e)de, is given by [72] 


1 


Tint 


C 

2 



rc(e) 

7 2 e 2 



(Bl) 


where cr(e') is the photo-nuclear interaction cross section of a nucleus of mass Am p by a photon of energy e' in the 
rest frame of the nucleus. 

Detailed tables of cr(e') are available in CRPropa [83, 84]. We use the numerical tools provided at [47] to calculate 
the interaction times for the photon field given by Eqs. (Al) and (A2). 

For illustrative purposes, the cross section can be approximated by a single pole in the narrow-width approximation, 


<7{s ') = 7 r < 7 r , 


S{e' - 4es) > 


(B2) 


where cr res is the resonance peak, r res its width, and e' es the pole in the rest frame of the nucleus. The factor of 1/2 
is introduced to match the integral (i.e. total cross section) of the Breit-Wigner and the delta function [85]. 

The mean interaction time is obtained substituting Eq. (B2) into Eq. (Bl), 


C 77 77 res c res r res 


Tint(-E) 


4 7 2 J o % < £ ) ®( 2 7£"4b) 

r d£ ,, 

-a»(e). 

Je’ 2 -y £ 


c 77 cx res c res r res 


4 ^ 2 Je' M /2 7 

Substituting (Al) into (B3) yields: 

1 _ 1 | (E b /E)P +1 E^E b 

“ n { (1 - /3)/(l - a) [(E b /E) a+1 - ( E b /E ) 2 ] + (E b /E) 2 E > E b 

where 


Tfc = 


2 E b (1 - /3) 

C 77 C7 res A if Ip r res ’{!.{) 


and E b = 


4es A m P 

2en 


(B3) 


(B4) 


(B5) 


The parameters characterizing the photo-disintegration cross section are: a res « 1.45 x 10 ~ 27 cm 2 A, r res = 8 MeV, 
and e' es = 42.65A -0 21 (0.925A 2 - 433 ) MeV, for A > 4 (A < 4) [74], The parameters for the photo-pion production 
cross section are: cr res ~ 5.0 x 10 -28 cm 2 A, r res = 150 MeV, and c' res = (rn 2 A — m 2 )/(2m p ) ~ 340 MeV [ 66 ]. 


Appendix C: Propagation in the Source Environment 


In our simple model we consider interactions and escape of particles treating the source environment as a leaky 
box. If at a given time, to, N(t = t 0 ) = N 0 particles are injected at random into the source environment, then the 
number of particles N remaining in the source at any later time t changes as 


dN 

dt 


-—N- —N, 

7"esc ^"int 


(Cl) 
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where T esc and Tint are the escape and interaction times respectively. Integration yields the time evolution of N as 


N{t) = N 0 e- 


where 


TescTint 

T = - 

'Tesc T ^"int 


The total number of escaping particles is given by 


iv fisr . = 


r 00 i 

- N(t) dt = N 0 

Jto T esc 


'Tint 


'Tq sc T ^int 


= No 


1 ^"esc/’ 7 "int 


= N 0 fe: 


and likewise the number of particles suffering interactions is 


(C2) 


(C3) 


(C4) 


iVint = No Te ° c = No - 1 -S No /int , (C5) 

^esc "f +nt 1 T ^"int/ 'esc: 

with N esc + iVi n t = Nq. As can be seen, 7V esc and A r j nt depend only on the ratio of the escape and interaction times, 
but not on the absolute value of either of them. 

In the following we consider sources at steady state, i.e. sources which are active long enough to justify integrating 
to infinity in Eq. (C4) and for which the injected flux equals the escaping flux. Interacting particles constitute the 
source for secondary particles of lower mass number. 

Since the particle trajectory in the source is treated as a random walk starting from a random position, the escape 
time of a secondary does not depend on the time it was produced. Therefore we can apply Eqs. (Cl) to (C5) also to 
the secondary particle production, which greatly simplifies the equations with respect to previous analytic approaches 
that had been developed for the extra-galactic propagation of cosmic-ray nuclei [79-82]. 


1. Single-Nucleon Emission 


The basic principle of the analytic calculation can be best illustrated by firstly describing the case where interactions 
with the photon field lead to the knock-out of a single nucleon, 


A + 7 —» (A — 1 ) + n/p, 


(C 6 ) 


and the nucleon carries away a fraction of 1/A of the initial energy of the nucleus. This approach has been successfully 
applied to the photo-disintegration (PD) during the extra-galactic propagation of nuclei (see e.g. [79]). It can also serve 
as a good approximation for the losses due to photo-pion production (PP) if nuclei are treated as the superposition of 
A individual nucleons (see e.g. [83, 84]). The interaction time is therefore the combination of the two processes, i.e. 


Tint — 


r PD 

'int 


+ T: 


PP 


r PD T PP 
'int 'int 


(C7) 


In this simplified propagation scheme, secondaries with mass A and energy E* originate from nuclei with energy 
E' = NEk E* and mass A + 1. They are produced at a rate 


Q(E*, A) = Q int 


A+l 


E*, A + 1 


dE’ 

dE* 


= Q 


A + l 
A 


E* , A + 1 I 7/int 


A + l 


E* , A+l 


A + l 


(C 8 ) 


where the factor 


dE' 


dE* 


is the Jacobian determinant needed to transform the differential injection rate from the primary 


to secondary energy. In analogy to Eq. (C4), a fraction of the secondaries escapes the source environment, 


g esc (£*, A) = Q(E*, A) Veac (E*, A), 


(C9) 


and the remaining particles interact again at a rate of 


Qint(E*, A) = Q(E*, A) 77 int (£*, A). 


(CIO) 
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This assumes that the escape probability of a secondary is independent of the time or position it got produced in 
the source environment. This calculation can be iterated to obtain the escape rate of any remnant with mass A* 
produced during the propagation of a nucleus of mass A!\ 


Q^(E*,A*,A') = Q(^E*,A' 


A' 

A* 


A' 

;(E*,A*) [] Tfint 

A°=A*+1 


A* 


E*, A < 


(Cll) 


The rate of nucleons being knocked out of nuclei during propagation via either of the considered processes i = PD/PP 
is 


Ql sc (E*,n+p,A') = Q 


A' E* A! „ / A+ E* 


r S /• 


A* =2 


A! 

n riu,t 

A°=A + 


/ A° E* 

V K i 


,A< 


(C12) 


with elasticities of the knock out nucleon given by kpd = 1 and Kpp = 0.8 and the fractional contribution from PD 
and PP given by 


/pd = 


1 


-) and /pp = 1 — /pd = 


1 


(C13) 


1 + rpo/rpp’ 1 + tpp/tpd 

The total escape rate of particles of mass A* from injected nuclei of mass A' is 

Q^t(E*, A*, A') = + + , (C14) 

where 5ai is the Kronecker delta. 


2. Branching Ratios from Photo-disintegration 


The propagation scheme described in the last section can be easily extended to take into account the emission of 
several nucleons or light nuclei in photo-nuclear reactions. We use the total interaction time and branching ratios for 
photo-dissociation from Talys [43, 44] as available in CRPropa and neglect multi-nucleon emission for photo-pion 
production since it can be safely neglected at the energies relevant here. 

Instead of the closed formulae derived above for the single-nucleon case, we now have the following recursive relation 
for the rate of produced remnant nuclei of mass A*. 


Q(E *, A*) 


A’-A* 

£ b(E u A* 1 A l )r lint {E l , Ai) Q(E U At) 

i=i 


dEj 

dE* 


(CIS) 


with Ai = A* + i, Ei = At/A* E* and \ ^§k\ = Ai/A*. b(Ei , A*, A/) is the branching ratio that gives the probability 
that a nucleus of mass Ai with energy Ei will have a remnant mass of A* after the interaction. The n knocked out 
nucleons and nuclei are calculated the same way but replacing the branching fraction by nA*,Ai ), i.e. the 

probability to produce n fragments of mass A* , and summing over n. The rest of the calculation proceeds as in the 
case of single-nucleon emission. 


3. Proton Interactions 

Once nucleons are generated in photo-disintegration they are assumed to either escape immediately in the case of 
neutrons, or to interact further via photo-pion production. The average elasticity of this process is Kpp = 0.8 and 
corresponds to a shift in energy of Algi! = lgKpp as 0.1. Since we perform the calculation in logarithmic bins of 
this width, proton interactions can be treated similarly to photo-disintegration as a “trickle-down” of particles fluxes 
subsequently shifted by one energy bin. For the reaction p + 7 — » tt + + n , the neutron escapes and the interaction 
chain is finished. In case of p + 7 —» 7 r° + p, the secondary proton has a reduced energy; it may also interact again. 
The neutron, proton and positive pion fluxes in an energy bin k are calculated from the recursive relations 

Q(E* k ,n)' = Q(E * k , n) + (1 - b pp ) Vint (E* k+1 , p) Q(E* k+1 ,p )'/« (C16) 


Q(E* k ,p)' = Q(E*,p) + b pp 7?i„t (E* +1 , p) Q(E k+1 ,p)'/n 


(C17) 
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and 


Q(E* k ,ir+y = (1 - b pp ) mnt(E* +7 , p) Q(E* +7 ,p)'/(1 - tv) , (C18) 

where b pp 0.5 is the branching fraction of the process p + 7 —» 7 r° + p and the un-primed fluxes are the sum of the 
knocked-out nucleons from Eq. (C15) and primary protons. The offset of 7 in the equation for the pion flux is due to 
the energy shift of the pions, A\gE = lg(l — k pp ) 0.7. 


Appendix D: Cosmic Ray Production and Propagation in an Expanding Universe 


To compare the spectra obtained in the last section, the particles need to be propagated to Earth. The number 
of cosmic rays per unit volume and energy in the present universe is equal to the number of particles accumulated 
during the entire history of the universe and is comprised of both primary particles emitted by the sources and 
secondaries produced in the photo-disintegration process. Herein, the variable t characterizes a particular age of the 
universe and t h indicates its present age. We adopt the usual concordance cosmology of a flat universe dominated 
by a cosmological constant, with Ha ks 0.69 and a cold dark matter plus baryon component fi m ss 0.31 [ 86 ]. The 
Hubble parameter as a function of redshift z is given by H 2 (z) = 77o[H m (l + z) 3 + Ha], normalized to its value 
today, Hq = 100 ft km s -1 Mpc^ 1 , with ft ~ 0.68 [ 86 ]. The dependence of the cosmological time with redshift can be 
expressed via dz = —dt(l + z)H(z). The co-moving space density of cosmic rays ncR of mass A from a population of 
uniformly distributed sources with (possibly age-dependent) emission rate per volume Q{E ', A',t) is given by 


ncR(E,A,A') 


dNa r 
dEdV 


n t 

) 


H d&> AA jE',E,t) 

dE 


Q(E',A',t) £(i) dE 1 dt , 


(Dl) 


where dS 6 aa' I dE is the expectation value for the number of nuclei of mass A in the energy interval ( E , E + dE ) which 
derive from a parent of mass A! and energy E' emitted at time t, and £(f) is the ratio of the product of co-moving 
source density and Q(E r , A! , t), relative to the value of that product today. Note that d£Z AA '/dE includes propagation 
effects both at the source environment and en route to Earth. 

We assume that the emission rate of cosmic rays is the same for all sources and the spectrum and composition is 
independent of the age of the universe, so that evolution of the volumetric emission rate with cosmological time can 
be described by an overall source evolution factor, £(t) discussed below. (It need not be specified whether this is due 
to an evolution of the number of sources or their intrinsic power.) We further assume, as per usual practice, that 
emission rate is fairly well described by a power-law spectrum. Under these general assumptions the source emission 
rate per volume takes the form 


C (^ eo (g) T e xp (-^-), (D2) 

where E ^ ax is the maximal energy of emitted protons, i.e., maximum rigidity of the accelerator, Z' is the nucleus’ 
atomic number, E 0 is some reference energy, and 


Qo 


n 0 


dN., 

~d& 


n 0 


dN A , 

dE' dt 


E'=Eq 
E' = Eq 


for bursting sources 
for steady sources 


(D3) 


Here, no is the number of bursts per unit volume per unit time and dN A i/dE' is the spectrum of particles produced 
by each burst, or for a steady source no is the number density of sources at 2 = 0, and dN A '/dE'dt is the UHECR 
production rate per unit energy per source. The cosmic ray power density above a certain energy E' min is given by 



(D4) 
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where T denotes the upper incomplete gamma function. 

The cosmological evolution of the source density per co-moving volume is parametrized as 


n s {z ) = n 0 £0) 


(D5) 


with £(2 = 0) = 1. We adopt for the fiducial model that the evolution of sources follows the star formation rate with 




(1 + z) a 
! + [(! + z)/b\ c 


(D6) 


where a = 3.26 ± 0.21, b = 2.59 ± 0.14 and c = 5.68 ± 0.19 [60]. Additionally we consider the family of evolution 
models parameterized as £(z) = (1 + z) m . 

To propagate the particles escaping the source environment to Earth we use the CRPropa framework [83, 84], 
For this purpose, we generate a library of propagated nuclei with A* = 1... Aj„ ax injected uniformly in light-travel 
distance. The latter corresponds to a non-evolving source distribution in comoving distance after accounting for the 
cosmological time dilation. We simulated particles up to Aj^ax = 56. Given this library of simulated particles, we can 
construct the propagation matrix for arbitrary source evolutions for each nuclear mass A * escaping the source 

and secondary mass A v at Earth. The elements of the propagation matrix give the expected number of secondaries in 
an energy interval [lg Ej , lg Ej + A] at Earth originating from nuclei at the source at an energy [lg E *, lg E* + A] for 
a given source evolution £(z(f)) and a uniform logarithmic spacing in energy with A = 0.1. Numerically, the elements 
are constructed via discretization of Eq. (Dl) 


A' 


AV 


n C K{E j ,A v ,A!) = 2 2 2 aE- 

A*=A„ i=j o=0 J 


wva Q^l(Et,A;,A') ata) a t a A Ef, 


(D7) 


where At a = tn/na and 


AV, 


ijliva 


1 Njffi%(E*,E* + A E*-E j ,E j + AEf, A* ■ A v \t a ,t a + A t a ) 


AR¬ 


AB,- 


N?™(E*, E* + A E *; A*; t a ,t a + At a ) 


(D8) 


For a non-evolving injection rate per unit volume, the number of generated events per bin is constant, Nf™ = Kf° n . 
Then, for any source evolution £[x:(i))], (D7) can be rewritten as 


A' 


\ p* "a 7yEarth 

n C n(E J ,A v ,A')= 2 2 A ') 2 iTW At a 
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A E, 
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. . A Ej 

A*=A„i=3 J 


(D9) 




where Yi P =’o" £[z(R>)] denotes the ([-weighted sum over all events generated with ( A*,E *) arriving at Earth with 
(A„, Ej) and Nf° n is the total number of generated events with (A*, Ei). Note that if binned in Azb = z max /b , then 
Nfub\Azb/Atb\ = constant, and hence (D9) can be rewritten as 


A' 


n C R.{E u A v ,^) = 2 


,A E* Yr a l 0 N^{z b ) 


A* =A U i = 3 TTaloKa 


AZq 

At a 


(DIO) 
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where \Az b /At b \ = (1 + z b )H(z b ) and z max = A z b n b . 

For a given spectrum of injected nuclei of mass A! we obtain the space density of cosmic rays at Earth with energy 
E and mass A, 


ncR.(E,A,A') = -^5_. 
v ; dEdV 


(Dll) 


For an isotropic arrival direction distribution (which is an excellent approximation based on current observations) the 
relation between the spectrum and the cosmic ray density is 




The total flux at earth of particles of mass A v is 


A' 


J{E,A V ) = ^ fiA’JJ&A^A’J, 

A'. = A„ 

where f(A') denotes the fraction of particles of mass A' injected at the source. 


(D12) 


(D13) 


Appendix E: Neutrino and Photon Production 


The results of the last two sections can be readily applied to obtain the flux of neutrinos at Earth from the decay 
of neutrons and charged pions. We approximate the emission rate of pions from photo-pion production by using 
Ktt = 1 — Kpp in Eq. (C12). The energies of neutrinos escaping from the source are given by the kinematics of the 
two-body decay of pions and the subsequent muon decay which we treat approximately by assigning a third of the 
muon energy to each of the decay products. In this way we construct a propagation matrix for pions, with 

j, = (y^, v e [pe)) and ^ = 7T“. Similarly, a propagation matrix for neutrons is obtained with v = (P e ) and M = n. 

Neutrino oscillation over astronomical distances modifies the initial flavor distribution of fluxes, in 

calculable ways. The relevant parameters for such a calculation are the three Euler rotations ( 612 , $23 > #13) and the 
CP-violating Dirac phase S. The current best-fit values as well as the allowed ranges of the mixing parameters at 
the la level are: 0 12 /° = 33.57±g;^, 6» 23 /° = 41.9±g;| © 50.3±|;|, 0 13 /° = 8.73+°;||, 6 /° = 2661“ [87]. The mixing 
probabilities are given by 


Pu 


— c 13 s 23 + ( c 12 c 23 + s 12 s 13 s 23 — ^ c 12C23Sl2Sl3S23Cs) 2 + (c 23 Si 2 + Ci 2 S 33 S 23 + 2C12C23S12S13S23C5) 2 > (El) 


Pv e ++Vp — 2 c 13 {c 12 S 12 c 23 T (c 12 + S 12 ) s 13 s 23 + C12 Si2 C23 S23 C$ (c 12 S 12 ) Si 3 } , 


(E2) 


Pve++v T = Pv e ++ 1^(023 $23 + 7t/2) , Pv T ^>v T = Pv^v^ifi 23 -* $23 + 7 t / 2 ) , (E3) 

and the unitarity relations 

P = 1 p - p p = l — P p p = 1 _ p — p (E4I 

where Cij = cos 0,j , Sij = sin 6 l:j , and cs = cos <5 [88]. The measurable neutrino flux at Earth is given by 

/ $ e \ / 0.55 0.24 0.21 W \ 

= 0.24 0.37 0.38 . (E5) 

\<S>t J \ 0.21 0.38 0.41 J \ J 

In addition to neutrinos, photons are produced from 7r° production and decay [89], and by photo-disintegration 
of high-energy nuclei followed by immediate photo-enrission from the excited daughter nuclei [90]. The 7-rays, elec¬ 
trons, and positrons produced in the decay of 7r° and it- trigger an electromagnetic (EM) cascade on the cosmic 
microwave background, which develops via repeated e + e~ pair production and inverse Compton scattering. Other 
contributions to the cascade are provided by Bethe-Heitler production of e + e~ pairs and 7-rays emitted during the 
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photo-disintegration process, after the photo-dissociated nuclear fragments de-excite. The net result is a pile up of 
7-rays at GeV < E 1 < TeV, just below the threshold for further pair production on the diffuse optical backgrounds. 
The EM energy then gets recycled into the so-called Fermi-LAT region, which is bounded by observation [91, 92] 
to not exceed w cas ~ 5.8 x 10~' eV/cm 3 [93]. The latest Fermi-LAT lirnts [92] do not significantly influence the 
determination of the w cas upper bound of [93] , because that bound is not very sensitive to the high energy bins added 
by [92], 

The photons coming from photo-pion production in the source environment were shown to be below the Fermi-LAT 
bound in Section III A. To place a bound on the contribution of photons from nuclear de-excitation to the Fermi- 
LAT diffuse gammas, without performing an explicit calculation, we can turn to the estimate of [57] which found 
uj cas aa 1.1 x 10“' eV/cm 3 assuming the high-energy IceCube spectrum to be entirely due to neutron beta-decay. Since 
in our model the neutrino flux from neutron decay is significantly below the IceCube spectrum, the corresponding 
de-excitation photon contribution to the Fermi-LAT data must be far below the limit. 
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